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ABSTRACT 

We perform a series of cosmological simulations using Enzo, an Eulerian adaptive- mesh refinement, 
N-body + hydrodynamical code, applied to study the warm/hot intergalactic medium. The WHIM 
may be an important component of the baryons missing observationally at low redshift. We investi- 
gate the dependence of the global star formation rate and mass fraction in various baryonic phases 
on spatial resolution and methods of incorporating stellar feedback. Although both resolution and 
feedback significantly affect the total mass in the WHIM, all of our simulations find that the WHIM 
fraction peaks at z ~ 0.5, declining to 35-40% at z = 0. We construct samples of synthetic O VI 
absorption lines from our highest-resolution simulations, using several models of oxygen ionization 
balance. Models that include both collisional ionization and photoionization provide excellent fits to 
the observed number density of absorbers per unit redshift over the full range of column densities 
(10^^ cm~^ < Nqyi ^ 10^^ cm~^). Models that include only collisional ionization provide better fits 
for high column density absorbers (A/'ovi ^ 10^^ cm~^). The distribution of O VI in density and 
temperature exhibits two populations: one at T ~ 10^-^ K (collisionally ionized, 55% of total O VI) 
and one at T ~ 10^"^ K (photoionized, 37%) with the remainder located in dense gas near galaxies. 
While not a perfect tracer of hot gas, O VI provides an important tool for a WHIM baryon census. 
Subject headings: cosmology: observations — intergalactic medium — quasars: absorption lines 



1. INTRODUCTION 

It is well established that the fraction of baryons in 
the universe that are easily observable drops from nearly 
100% at high redshift to less than half by the cur- 
rent epoch. Predictions from Big Bang Nucleosynthe- 
sis (Olive et al. 2000; Steigman 2007) and measurements 
of the cosmic microwave background by Komatsu et al. 
(2010) are in agreement with the cosmic baryon budget 



measured by observations of the Lya forest at 2: > 4. Yet 
baryon surveys at z = consistently finds that as many 
as 60% of those baryons are "missing" (Per sic fc Salucci 
1992; Bristow & Phi llipps||1994| ^^Fukugita et al. 1998; 
Fukugita & Peebles 2004|) . In the low-redshift intergalac- 
tic medium (IGM), approximately 30% of the baryons re- 
side in photoionized diffuse Lya absorbers (the "Lyman- 
alpha forest") and 10% in hot gas traced by O VI 
absorbers (Danforth & Shull 2008). For reviews of the 
missing-baryons problem and observational attempts at 
a solution, see Bregman (2007) and Shull (2003). Early 
numerical s imulations ( |Cen fc Qstriker 1999^ Dave et alT] 



1999[ 2001 ) predicted that a majority of the missing 
baryons are located in a gaseous phase of moderate to 
low density and temperatures lO^-lO*" K, known as the 
warm/hot intergalactic medium (WHIM). These simula- 
tions and many that followed showed that many baryons 
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are heated into the WHIM phase via gravitational shocks 
from gas falling onto dark-matter filaments and halos, 
and through galactic outfiows and stellar feedback. 

WHIM gas is quite tenuous and difficult to observe, 
with baryon overdensities of typically 5h = pb/pb ^ 
0.1 — 100. Here, pb and pb are the baryon density and 
mean baryon density, estimated at p^ = ftb Pcr,o (1 + 
z)^ = (4.26 X 10"^^ g cm"^)(l + z)^, where pcr,o is the 
critical density of the universe at z = 0. The method 
for probing the WHIM that has proved most fruitful 
thus far is the detection of the O VI doublet (AA1032, 
1038) in absorpt ion in the low -redshift Lya forest, first 
accomplished by Tripp et"aL] ([2000 ) . In collisional ion- 
ization equilibrium (CIE), O VI reaches its pea k abun- 
dance fraction, f ovi ~ 0.2, at T ~ 10^'^^ K (Sutherland 



& Dopita 1993), potentially making it an ideal tracer 
of WHIM gas below 10^ K. Since the first detections, 
large surveys of low-redshift O V I absorption lines have 
been completed ([Da nforth fc Shuh 2005', '2008'; 'Danfo rth| 



et al. 



2006 Tripp et al. 2008; Thom & Chen 2008a| 



providin g statistical datasets for comparison with simu- 



lations. Furlanetto et al. (2005) showed that O VI and 
O VII absorption statistics can be fit reasonably well 
by models that assume these systems to be generated 
in a network of virial and infall shocks that surround 
cosmological halos. Cen fc Fang (2006) were able to 
broadly reproduce the observed distribution in column 
density of O VI absorbers, quantified as the number 
of absorbers per unit redshift (dJV/dz)^ by including a 
non-equilibrium treatment of oxygen ionization balance 
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and stellar feedback from galactic superwinds. However, 
they did not explicitly study whet her these absorbers are 
tracing the WHIM. More recently, |Oppenheimer fc Dave 
f2009) found that a majority of O VI absorption systems 
in their simulations arise from cold (T ~ 15, 000 K) gas 
that is primarily photoionized. They suggeste d that O VI 
may n ot trace the shock-heated WHIM at all. Kang et al. 
( 2005 ) arrived at a similar conclusion. However, t he raT 
diative cooling method of [Oppenheimer fc Dave| ( 2009 ) 
only accounted for the intluence ot photoionization on 
the cooling of H and He, but not the metals, whose cool- 
ing rate was calculated assuming CIE. At moderately 
high metallicities {Z > O.IZ©), this can lead to an over- 
estimate of the total cooling by more than an order of 
magnitude for T < 10^ K. An exce llent demonstration o f 
this is given in Figure 10 of Tepper-Garcia et al.| (|2010|). 

In this work, we present the results of a new s^t of nu- 
merical simulations designed to study the nature of the 
WHIM and its relation to low-redshift metal absorption 
systems. This paper is the first in a series in which we use 
the results of these and additional simulations to study 
the observational signatures of the WHIM. Our primary 
goals in this first paper are to validate our approach and 
provide initial predictions. In Section [2j we detail the 
numerical methods employed in this work. In Section 
|3j we discuss the initial results of these simulations. To 
understand the extent to which these simulations can be 
trusted, we begin by examining how well they are con- 
verged, a topic that has received little previous attention. 
We use our highest-resolution simulations to create sam- 
ples of synthetic O VI absorption lines for comparison 
with observations. We then take a closer look at the 
physical environment associated with the O VI in our 
simulations, in order to gauge its usefulness as a WHIM 
tracer. We conclude with a discussion and summary in 
Section [H 

2. METHODS 

In order to understand the dependence of WHIM char- 
acteristics on the properties of the simulations, we carry 
out a large number of cosmological simulations with vary- 
ing spatial resolution, box size, and physical parameters. 
In every simulation, the number of grid cells is equal to 
the number of dark-matter particles. All of the simula- 
tions are initialized at z = 99 in a ACDM universe with 
a power spe ctrum of density fluctuations given by Eisen- 
stein fc Hu (1999). The cosmological parameters are (1^5, 
i^cDM, ^lA.ti^, n) = (0.0441, 0.2239, 0.732, 0.704, 0.82, 
1), in reasonable agree ment with the WMAP-7 findings 
(Komatsu et al. 2010). All simulations with the same 
box size are mitialized using identical random seeds, and 
thus have the same large-scale structure, allowing us to 
isolate the effects of resolution. Simulations with differ- 
ent box sizes have different random seeds. A summary of 
the simulations and their distinguishing features is given 
in Table [l] We now describe the setup and methodology 
in further detail. 

2.1. The Enzo Code 

We performed the entire suite of simulations with 
the Eulerian adaptive mesh-refin ement (AMR), hydrody- 



namics + N-body code, Enzcj^ ( [Bryan fc Norman 



1997 



O'Shea et aH2005a|b ). The equations of hydrodynamics 
are solved u s ing th e iZEUS hydrodynamic solver of |Stone| 
fc Norman (11992*). The evolution of collisionless dark 



matter is computed with an N-bod y solver that uses a 



particle-mesh ( PM) gravity solver (Efstathiou et al.|1985 



Hockney fc Eastwood 1988). The AMR functionality o: 



Enzo provides the capability to dynamically add reso- 
lution to regions of the computational domain where it 
may be required for a variety of reasons, such as baryonic 
or dark-matter overdensity, resolving the Jeans length or 
shocks. This makes it possible to have high resolution 
in regions of interest while ignoring less important areas, 
thus limiting computational costs. However, since the 
WHIM is thought to reside in regions of relatively low 
overdensity (0.1 < Sh = Pb/ Pb ^ 100)? a very large frac- 
tion of the computational domain must be properly re- 
solved. Rather than allowing for vast regions of adaptive 
refinement throughout the domain, which can introduce 
significant computational cost simply to manage to the 
hierarchy of grid patches, we instead choose to run with 
AMR disabled and to perform simulations with large, 
static grids, referred to here as "unigrid" simulations. 

2.2. Star Formation and Feedback 

Star formation is a process that occurs on physical 
scales many orders of magnitude below the resolution 
limit of cosmological simulations, yet the products (ra- 
diation, thermal energy, and heavy elements) have sig- 
nificant influence on the evolution of structure on large 
scales. Numerical simulations must employ models based 
on our limited understanding of the physical conditions 
that produce stars and the impact they have on their 
surroundings, on the physical scales relevant to the sim- 
ulations. Numerous prescriptions for feedback have been 
introduced with s lightly varying parameterizations, (e.g . , 



Cen fc Ostrikerl 



Schaye & Dalla 



j ligntly varying pai 
([1992I; [Sprh^ge H 

Vecchia (2008)), 



Hernquist (2003a 
that are capable of 



broadly reproducing the observed star formation history 
[opkins fc Beaconi||2006 ) . 

We use a modi tiea version of the prescription of |Cen| 



Ostriker (1992) in our simulations, which we describe 
briefly here. Any grid cell is capable of forming stars 
if the following criteria are met: baryon overdensity, de- 
fined here and only here as pb/{pc,o{^^z)^) > 100, where 
Pc,o is the critical density at 2; = 0; the velocity diver- 
gence is negative; and the cooling time is less than the 
dynamical time. In all other places, we define baryon 
overdens i ty as ph/pb- The original prescription of Cen fc 
Ostriker (1992) calls for the additional check of Jeans in- 
stability. However, we find that in large-scale, fixed-grid 
simulations, this final criterion is always met, following 
satisfaction of the first three requirements. Therefore, 
we ignore the check for Jeans instability as a cost-saving 
measure. If all three listed criteria are satisfied, then a 
"star particle," representing a large collection of stars, is 
formed within the grid cell with a total mass. 



At 



(1) 



yn 



http://code.google.eom/p/enzo 



where is an efficiency parameter, mceii is the baryon 
mass in the cell, t^yn is the dynamical time of the com- 
bined baryon and dark matter fluid, and At is the hydro- 
dynamical timestep. Subsequently, this much gas mass 
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Table 1 

Parameters of Simulations 



Run 


(Mpc//i) 


1 /Q 

^^cells 


Ax 
(kpc/h) 


rridm 
(Mo/Zi) 


Feedback 
Method 


Metal 
Yield (y) 


OK ot^a n 
ZD-ZdoJj 


OK 
ZO 


O K^? 

zOD 


OQ 


D X iU 


U 




OK QQ/l n 

zO_oo4_U 


OK 
ZO 


QQ/l 

oo4 


a K 
DO 


z X iU 


U 




OK Ki o n 


OK 
ZO 


K1 O 

Oiz 


A n 

4y 


< X iU 


U 




OK 'Ti^Q n 


OK 

zO 


(Do 


Q Q 
OO 


z X iU 


U 




OK 1 no/1 n 


OK 

zO 


1 no /I 
iUz4 


O/l 

z4 


y X iU 


U 




OK OK<5 1 


OK 

zO 


O K<5 

zOb 


OQ 


X iU 


i 


n noK 
U.UzO 


OK QQ/l 1 


OK 

zO 


QQ/l 

oo4 


<5 K 
DO 


z X iU 


i 


n noK 
U.UzO 


OK K1 O 1 


OK 

zO 


K 1 O 

Oiz 


A n 

4y 


/ X iU 


i 


n noK 
U.UzO 


OK 'Ti^Q 1 


OK 

zO 


(DO 


Q Q 
OO 


z X iU 


i 


n noK 
U.UzO 


OK 'Ti^Q 1 K 


OK 

zO 


(Do 


Q Q 
OO 


z X iU 


i 


n noK 
U.UUO 


OK OK<^ O 
ZO_ZOO_Z 


OK 

zO 


O K<^ 

zOO 


HQ 

y© 


D X iU 


o 
z 


n noK 
U.UzO 


OK QQ/l O 


OK 

zO 


QQ/l 

oo4 


<5 K 

Do 


z X iU 


o 
z 


n noK 
U.Uzo 


OK K1 O O 

zO_Oiz_z 


OK 

zO 


K 1 O 

Oiz 


/I n 

4y 


/ X iU 


o 
z 


n noK 
U.Uzo 


OK TfiQ O 


OK 

zO 


(Do 


Q Q 
OO 


z X iU 


o 
z 


n noK 
U.UzO 


Kn K 1 o n 


50 


512 


98 


6 X 10^ 







50_768_0 


50 


768 


65 


2 X 10'^ 







50_1024_0 


50 


1024 


49 


7 X 10^ 







50_512_1 


50 


512 


98 


6 X 10'^ 


1 


0.025 


50_768_1 


50 


768 


65 


2 X 10'^ 


1 


0.025 


50_768_lb 


50 


768 


65 


2 X 10'^ 


1 


0.005 


50_1024_1 


50 


1024 


49 


7 X 10^ 


1 


0.025 


50_512_2 


50 


512 


98 


6 X 10'^ 


2 


0.025 


50_768_2 


50 


768 


65 


2 X 10'^ 


2 


0.025 


50_768_2b 


50 


768 


65 


2 X 10'^ 


2 


0.005 


50_1024_2 


50 


1024 


49 


7 X 10^ 


2 


0.025 



1 /3 

Note. — Run: simulation name, given as (O-(Keils)-(feedback method), 

1 /3 

where ^ is the comoving box size and N^^^^^ is the number of grid cells along 
each edge of the box. The total number of grid cells is equal to the number 
of dark matter particles. Ax: comoving grid-cell size, rridm' dark matter 
particle mass. Feedback method: 0: adiabatic (no radiative cooling, feed- 
back from star formation, or ionizing UV background); 1: local method; 2: 
distributed method. Metal yield: fraction of stellar mass returned to gas as 
metals (mmetais/^*)- 



is also removed from the grid cell as the star particle is 
formed, ensuring mass conservation. 

While the star particle is formed instantaneously 
within the simulation, feedback is assumed to occur over 
a longer time scale, which more accurately reflects the 
gradual process of star formation. In a single timestep of 
the simulation, the amount of mass from a star particle 
converted into newly formed stars is given by 

Am,f = — e-(^-t*)/ta.n ^ (2) 

^dyn ^dyn 

where t is the current time and the formation time 
of the star particle. Stellar feedback is represented by 
the injection of thermal energy and the return of gas and 
metals to the grid, all in amounts proportional to Amgf. 
The thermal energy, e, and metal mass, mmetais^ returned 
to the grid are given by 

e Amsf e , mmetais Amsf y , (3) 

where c is the speed of light, e is the ratio of thermal 
energy output to total rest-mass energy of the star par- 
ticle, and y is the fraction of stellar mass returned to 
the grid as metals. As per Equation |2j the rate of star 
formation, and hence the rate at which thermal energy 
and metals are injected into the grid by the star particle, 
rises linearly, peaking after one dynamical time, then de- 



cays exponentially after that. For a given star particle, 
30%, 60%, and 90% of the total feedback production has 
occurred by 1, 2, and 4 dynamical times, respectively. 
For a threshold baryon overdensity (as defined above) of 
100, corresponding roughly to a total overdensity of 600 
(assuming the cosmic ratio of baryonic to dark matter), 

tdyn-l^^^ (l + ^)-3/2 Gyr. 

We choose parameters for the prescriptions of s tar for- 
mation and feedback similar to those adopted by |Cen fc| 
Ostriker (1992). We assume a star-formation efficiency, 
= 0.1, a 25% fraction of stellar mass returned to the 
grid as gas, and e = 10~^. Of the gas returned to the 
grid, 10% is returned in the form of metals, for a total 
metal yiel d, y = 0.025, consistent with the calculations 
of Madau et al.| |T996). This metal yield, y = 1/40, is 
consistent with average values in the Milky Way, with a 
mean star-formation rate (SFR) of ~ 3 Mq yr~^, a core- 
collapse supernova rate of 1 SN every 40 years (from pul- 
sar counts), and an IMF-averaged metal yield of ~ 3 Mq 
per supernova. For a Salpeter IMF, the ratio of super- 
nova en ergy to metal eject a rest-mass energy is ~ 0.0015 
( Bookbinder et al.| ^1980). With our choice of y and e, 
this implies that approximately 27% of the total super- 
nova energy is converted into feedback energy. As one of 
the main goals of this work is to understand the relation 
between metals and the WHIM, we have also run three 
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additional simulations with the metal yield decr eased by 
a factor of 5 to y = 0.005. Schaye et al. (2010) investi- 
gated the dependence of the global star formation history 
on the various star formation and feedback parameters, 
finding it difficult to induce a significant deviation. A pa- 
rameter study performed as a part of this project showed 
the star formation history to be robust to changes in the 
prescription parameters. This effect is primarily due to 
the fact that star formation is limited by the amount of 
cold, collapsing gas, which is independent of the choice 
of parameter values. 

In this paper, we investigate two different methods for 
injecting feedback into the grid: (1) depositing all of the 
feedback into the grid cell where the star particle exists; 
and (2) distributing an equal amount among the central 
cell and its 26 nearest neighbors. Since this entire proce- 
dure is designed to account for processes unresolvable by 
the simulation, neither method should be considered to 
be representative of the actual physical processes that are 
occurring, but are, rather, an attempt to mimic the end 
state of those processes at the multi-kpc scale. Injecting 
feedback into a single cell has well-known practical limi- 
tations: in particular, the great increase in temperature 
produced by the injection of kinetic or thermal energy in 
a grid cell with reasonably high density can lead to cool- 
ing times that are much shorter than the hydrodynami- 
cal timestep, and thus, over-cooling. This diminishes the 
ability of the stellar products to propagate away from the 
source. This issue can be mitigated by making the as- 
sumption that feedback is able to spread to neighboring 
grid cells before energy is lost through radiative cooling. 
This technique was employed by Cen fc Fang, ( ,2006 ) to 
mimic the effect of galactic superwinds. In this work, we 
distribute feedback to the neighboring cells in the form 
of thermal energy, gas, and metals, and we do not find 
it necessary to add kinetic energy. Hereafter, we refer to 
the method of injecting stellar feedback into the central 
cell and its 26 neighbors as the distributed method, and 
the method of inserting feedback solely into the central 
cell as the local method. 

2.3. Radiative Heating and Cooling 

Of utmost importance to properly modeling the evo- 
lution of the WHIM is an accurate determination of 
the thermal state of the baryons within the simulation. 
Wiersma et al. (2009) showed that including both the ra- 
diative cooling from metals and the pho toelectric heating 
from the UV metagalactic background (Haardt & Madau 
2001 ) is vital for simulations of large-scale structure. Un- 
fortunately, a direct calculation of the heating and cool- 
ing rates, including metals, requires solving a network 
of rate equations that is far too large to be viable for 
millions, and now billions, of grid cells. 
We use an efficient computationa l method, origina ted 



by Smith et al. (2008) and similar to 
coupling a multidimensional table o: 
values , pr ecomputed with the p hotoionization software , 



Shen et al.| ( |2010D of 

heating and cooling 



values , pr ecomputed witn tne p notoionization software , 
Cloud^r] (Ferland et al. '19981 , to a non-equilibrium, 
primoroial cheni istry network (Abel et al. |1997 Anm 



nos et al.| |1997[ ) . The primordial chemistry network 
solves for the evolution of six species (H, H+, He, 
He+, He++, and e~), including collisional ionization and 

^ http://nublado.org 



photo-excitation/ionization rates, using a backward dif- 
ferencing formula. We use Cloudy version 07.02.01 to 
construct a grid of heating and cooling rates that vary 
with density, metallicity, electron fraction, redshift, and 
temperature. The redshift determines the natur e of the 
UV m etagalactic background, given by Haardt & Madau 



( 2001[) . We construct an analogous grid of values run 
with H and He and subtract this from the larger grid, 
so that only the metal heating and cooling rates remain. 
The rates are normalized by the H number density times 
the electron number density {nn x Uq). During the sim- 
ulation, we use the electron density from the primordial 
network plus a small, metallicity-dependent correction 
to calculate the heating and cooling contributions of the 
metals. We turn on the ionizing (UV) background at 
z = 7. The primordial chemistry network follows the 
photoionization of H and He, and we interpolate from 
the grid of data made with Cloudy to include both the 
heating and cooling due to the metals. For z > 7, we use 
a similar grid of heating and cooling data, created with 
no radiation background, assuming collisional ionization 
only. The metal heating and cooling terms are included 
within a coupled chemistry and cooling solver. In order 
to maintain accuracy, the solver subcycles the cooling 
and chemistry routines together on timesteps that are 
no larger than 10% of formation/destruction timescales 
{pi/ pi) of e~ and H, 10% of the cooling time (e/e, where 
e and e are the internal energy and its time derivative), 
and half of the hydrodynamic timestep. Coupling of the 
chemistry and cooling solver is particularly crucial in re- 
gions of high density and metallicity, where the cooling 
time is very short and the gas can transition from ionized 
to neutral within a single hydrodynamic timestep. 

3. RESULTS 

3.1. The Star Formation History 

3.1.1. The Effects of Radiative Cooling 

We first attempt to quantify the effect of the radiative 
cooling method chosen by comparing it with a set of four 
simulations run with less complete treatments. Each of 
these test simulations has a 25 Mpc//i box with 512^ 
grid cells and dark matter particles, with stellar feed- 
back injected into only one cell. These simulations are 
compared to run 25_512_1. In Figure [l] we plot the co- 
moving star formation rate density for run 25_512_1 and 
all of the test simulations, along with the comp ilation of 
observational dat a on comoving SFR density by |Hopkins| 
Beacom (|2006). The results of this experiment can be 



summarized very simply: a higher cooling rate results in 
a higher SFR. The removal of radiative cooling from met- 
als results in a significant decrease in the SFR from z ~ 5 
to the present. In the absence of both metal cooling and 
the radiation background (purple line in Figure [IF, the 
peak in the star formation rate is lower by a factor of a 
few than the control run, and occurs slightly earlier (just 
before z = 2, as opposed to just after.) The low-redshift 
slope is also much steeper, resulting in a SFR at z = 
that i s lower by an order of magnitude. [Schaye et al. 
( 2010[) observe a similar reduction in the peak SFR, but 
they do not see the steepening of the slope at low redshift 
reported here. 

The UV background provides additional heating to the 
gas, making it more difficult to form stars. With only pri- 
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0.4 0.6 
log (1 + z) 

Figure 1. Effect of radiative cooling and radiation background 
on the comoving star formation rate density plotted over the ob- 
servational data compiled by Hopkins Sz Beacom ( 2006). Shown in 
black is run 25_512_1, with the tull coolmg model (non-equilibrium 
H/He cooling plus a radiation background and metal cooling that 
includes the effect of the radiation background. Also shown are 
four additional test simulations with 25 Mpc//i box size and 512^ 
grid cells with less complete cooling methods. Purple: only the 
cooling from atomic H/He with no radiation background; blue: 
cooling from H/He with no radiation background and metal cool- 
ing that does not include the effects of the radiation background; 
green: cooling from H/He only plus a radiation background; red: 
cooling from H, He plus radiation background and metal cooling 
that does not account for the radiation background. 

mordial cooling, the addition of the ionizing background 
has a far greater influence on the ability of the gas to 
form stars than when metals are included. In the simu- 
lation with the radiation background but no metal cool- 
ing (green line in Figure [T]), the SFR is more than half 
an order of magnitude lower than the simulation with no 
radiation background and no metal cooling, or over one 
and a half orders of magnitude below the control run, 
which has both a radiation background and metal cool- 
ing. When metal cooling is included, the presence of the 
radiation background makes much less of a difference, as 
can be seen with the blue and black lines in Figure [l] 
Photoionization is more effective at reducing the cooling 
efficiency of primordial gas because of the relative ease 
with which H and He (compared with heavier elements) 
are fully ionized, thus eliminating the ability to cool via 
line emission. Because primordial coolants alone are re- 
duced more efficiently by the presence of the radiation 
background, it is much more difficult for the gas in these 
simulations to continue to cool and form stars. Finally, 
if the metal cooling method does not account for the in- 
fiuence of the radiation bac kground, as was neglected in 
[Oppenheimer & Dave (2009|, the potential exists for a 
significant overestimate in tne total cooling. This may 
have been the primary reason that the above work re- 
ported the existence of much colder O VI absorbers than 
had previously been claimed or is seen in this work. The 
red line in Figure [l] represents a simulation in which the 
metal cooling rates do not reflect the presence of the ra- 
diation background. Interestingly, we do not see such a 
dramatic influence on the SFR, as is evidenced by the 
nearly identical red and black lines in Figure [l] 

Not all simulations treat the radiative cooling rate con- 
sistently with the ionizing background and metal trans- 
port. Because the peak of the radiative cooling rate, at 
temperatures 4.5 < logT < 6.0, arises primarily from 



electron-impact excitation of bound states of abundant 
metal ions (C III, C IV, O VI, Ne VIII), a strong UV 
radiation field will photoionize away those coolants and 
decrease the cooling. As can be seen in Figure Jl] our 
most complete cooling model over-predicts the SFR at 
z < 2 m this series of simulations, owing to the inabil- 
ity of the feedback prescription to drive the met als away 
from their source ( Springel fc Hernquist||2QQ3 a1). As dis- 
cussed below, we are able to bring the predicted SFR 
into better agreement with observations simply by dis- 
tributing the thermal energy and metals over multiple 
grid cells. 

3.1.2. Convergence 



'Schaye et al. (2010) note the difficulty of matching the 
observed ISI^'K at both low and high redshift within a 
single simulation. At high redshift, extreme resolution is 
required to resolve the low-mass halos that dominate the 
global SFR. At low redshift, a large box size is needed 
to capture the most massive objects forming from large- 
scale perturbations in the primordial density field that 
are only recently becoming nonlinear. In Figure [2j we 
plot the SFR versus redshift and the corresponding stel- 
lar density evolution for all simulations listed in Table [l] 
with metal yield, y = 0.025. We do not see convergence 
in the SFR at high redshift, as noted in previous st udies 



dSpringel fc Hernquist||2003b[ [Schaye et ar||201Q] ). In- 

5lutic 



creasing the spatial and gas/ dark matter mass resolut ion 
which do not vary independently in our simulations, al- 
lows more low-mass halos to be resolved at progressively 
higher redshifts, continually increasing the SFR. Con- 
vergence at z = is achieved in our lowest-resolution 
simulation when using the local feedback method. When 
using the distributed method, convergence at z = is 
achieved in our second-lowest-resolution simulation. As 
the resolution is increased, the redshift at wh ich the SFR 
is converged also increases, as was noted in "Springe 
Hernquist (2003b). With the local feedback method, 
^FTToTThe two highest resolution runs (25_512_1 vs. 
25_768_1 and 50_768_1 vs. 50_1024_1) agree to within 
10% for z < 1.5. For the distributed- feedback simula- 
tions, this occurs slightly later, at z < 1.25. Since the 
interval from z = 1 ^ represents a majority of the age 
of the universe, the resulting stellar densities at 2; = are 
reasonably converged for all but the most poorly resolved 
simulations. 

3.1.3. Feedback and Box Size 

All simulations that employ the local feedback method 
overpredict the SFR at low redshift. This is due to the 
over-cooling that occurs when feedback is injected only 
into the central cell, unphysically dri ving up the gas t em- 
perature a nd the cooling rate (Katz et al. pTOQG; Baloghj 



et al. 2001j). The local feedback simulations with equal 
resolution have nearly identical SFRs at high redshift. 
However, those with the larger box size show slightly 
lower peak rates and a steeper decline at low redshift. 
This is likely due to the ability of the large-box sim- 
ulations to form structures with higher virial temper- 
atures, in which much of the gas becomes too hot to 
form stars. The simulations that use the distributed 
method display the best agreement with observations. 
In comparison to their local feedback counterparts, the 
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Figure 2. Top panels: comoving star formation rate (SFR) density versus redshift for simulations with varying resolution and box size 
plotted over the observational data compiled by Hopkins &: Beacom (2006). Simulations with the same resolution are shown with the same 
color: red, Ax = 98 kpc/h; green, 64 kpc/^blue, 49 kpc//i; and black, 33 kpc/h. Note that spatial and mass resolution are not varied 
independently in our simulations (see Table [T]). Simulations with 25 Mpc/h box sizes are shown with solid lines, and those with 50 Mpc//i 
box sizes are shown with dashed line s. Bottom panels: st ellar density in relation to critical density corresponding to top panels plotted 
over the observational compilation of |Wilkins et al.| ( |2008| . Left panels denote simulations with local feedback, while right panels denote 
those with distributed feedback. 

normal metal yield. To locate halos within our simula- 
tions, we use a version of the HOP halo- findi ng algo- 
rith m (Eisenstein fc Hut||1998 ), parallelized by ,Skory et 



distributed-feedback simulations have peak SFRs lower 
by 35-40%. 

The slope of the SFR with redshift (at z < 2) is also 
noticeably steeper in the distributed feedback runs. How- 
ever, unlike the local feedback simulations, there are no 
discernable differences between the small-box and large- 
box simulations of equal resolution when using the dis- 
tributed method. Comparing the highest resolution sim- 
ulations of each feedback method, we find the use of 
local method results in the formation of approximately 
1.7 times more stars than the distributed met hod. In 
comparison with observational measurements by jWilkinsl 
[et al.^ (2008) of the stehar mass density, our simulations 
still appear to pr oduce too niany stars by a factor of a 
few. However, ^Wilkins et al. (|2QQ8^ point out that in- 
dependent observations of the SFFTand the stellar mass 
density are not consistent with each other; the observed 
SFR implies a stellar density 2-3 times higher than what 
is seen. The final conclusion is that spreading stellar 
feedback over 27 cells, rather than just 1, is largely able 
to overcome the classic over-cooling problem, yielding a 
much better fit to fit to observations. 

3.2. The Halo Mass Function 

In Figure |3) we plot the cumulative halo mass func- 
tion for all the simulations with distributed feedback and 



al. 



(20l0p. We show only the mass functions for simula- 
tions with distributed feedback. For the range of halos 
captured in our simulations, there are no notable differ- 
ences for varying baryonic phy sics. In c omparison with 
the analytical fit of [Warren et al.| ( [200 6 ), we find our 
simulations to provide a good match for halos with ap- 
proximately 400 particles or more, regardless of particle 
mass. This result is in goo d agreement with numerous, 
prior Enzo simulations, e.g., O'Shea et al. ^ 
et al.| ( 12010 ); Skillman et al. (2010). 



(2005b); 



Trenti 



3.3. The Evolution of Baryon Phases 
3.3.1. Convergence 

In Figures [4] and |5| we plot the evolution of the baryon 
fraction in various IGM phases as a function of redshift 
for the majority of our simulations. We group the plots 
in rows and columns, showing box size and feedback pa- 
rameter so as to focus on convergence. We divide the 
baryons into four phases: warm (T < 10^ K, pb/pb < 
1000); WHIM (10^ K < T < 10^ K, pb/pb < 1000); hot 
(T > 10'^ K, pb/pb < 1000); and condensed {pb/pb ^ 
1000). We also include the total mass in stars in the 
condensed phase. In practice, the gas in the condensed 
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Figure 3. The cumulative halo mass function for all of the simulations with distributed feedback and normal metal yield. Simulations 
with 25 Mpc//i boxes are shown in the le ft panel while simulations with 50 Mpc/h boxes are shown in the right panel. Plotted in orange is 
the analytical fit of Warr en et al.| ( [2006| >. The colors shown here, which indicate the dark matter particle mass within the simulation, are 
consistent with Figure [2] 
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Figure 4. Evolution of baryon mass fraction in the warm (T < 10^ K), WHIM (10^ K < T < lO'^ K), hot (T > lO'^ K), and condensed 
{Pb/Pb ^ 1000) phases for the 25 Mpc/h box simulations. The condensed phase includes both gas with baryon overdensity 6h > 1000 and 
the total mass in stars. For the warm, WHIM, and hot phases, only gas with 6h < 1000 is included, such that the total is unity. The top 
row shows the adiabatic simulations, the middle row shows the simulations with local feedback, and the bottom row shows the simulations 
with distributed feedback. 



phase is almost entirely warm. If we define gas with 
Pb/pb < 1000 as diffuse^ we see that the ratio of dif- 
fuse to condensed WHIM gas is generally around 20:1 at 
z = 0. The masses of diffuse and condensed hot gas are 
within a factor of a few of each other, but neither con- 
tributes significantly to the total baryon budget within 
our simulations, mainly due to the somewhat small box 
sizes. Note that these phase cuts are slight l y diffe rent 
from those of other wo rk, such as jDave et al.| (|2001p and 



Cen & Ostriker ( 2006|) . However, the results are insensi- 
tive to these minor ditferences in definition. We have also 
measured the evolution of the baryon fractions, with the 



overdensity threshold dividing the condensed phase from 
the other three phases set to 200, and find the variations 
to be marginal. 

In Figure |4j we focus on the simulations with the 25 
Mpc/h box size. The top row of Figure |4j which shows 
the evolution of the adiabatic simulations, reaffirms the 
idea that shock heating from structure formation plays 
the primary role i n the evolution o f the WHIM ( |Dav6| 
^~al][T999| [2QQI| |Cen fc Ostriker] [l 999 J . Without ra- 
diative cooling, sta!r formation, and the LJV background, 
the general trends are similar. Warm, photoionized gas 
makes up roughly 90% of the total baryon mass at z = 4, 



8 



Smith et al. 



1.0 
0.8 
0.6 
0.4 k—- 
0.2 



Warm 



WHIM 



0.0 
1.0 
0.8 
0.6 

0.2 
0.0 
1.0 
0.8 
0.6 
0.4 
0.2 
0.0 



Adi.abati.c 



Local . 



Distributed 







1 2 3 4 
z 




0.03 
0.02 
0.01 



Hot 



Condensed 



- V"'-. 


Ax / h 

98 kpc _ 

- ■ - 65 kpc 
49 kpc 









Figure 5. Evolution of baryon fraction in the warm, WHIM, hot, and condensed phases for the 50 Mpc/h box simulations. The top row 
shows the simulations with local feedback and the bottom row shows the simulations with distributed feedback. Phase divisions are the 
same as in Figure [4] 



decreasing to ~40% by 2: = 0. This is offset primarily 
by an increase in WHIM gas from ~10% to 40-50% over 
the same interval. In the adiabatic simulations, the mass 
fractions of warm and hot gas evolve identically, inde- 
pendent of resolution. Since the box size in these simu- 
lations is rather small, massive halos with virial temper- 
atures Tvir > 10^ K are practically non-existent. Thus, 
the fraction of hot gas remains neglible at all times. The 
spike in the hot gas fraction at low redshift seen in all 
of the simulations shown in Figure [4] occurs due to a 
merger. One resolution effect that exists in the adiabatic 
simulations is a tradeoff between decreasing WHIM gas 
and increasing condensed gas with increasing resolution. 
As the resolution increases, the increased gravitational 
force resolution allows more low-mass halos to collapse 
and cross the critical overdensity threshold to be con- 
sidered condensed gas {pb/pb > 1000) before becoming 
pressure supported. In addition, higher resolution allows 
for smaller-scale perturbations to be resolved in the ini- 
tial density field, which will also increase the mass in 
collapsed objects. 

The addition of the full complement of physical pro- 
cesses to the simulations significantly alters the conver- 
gence properties. For both feedback methods, an in- 
crease in resolution yields a slightly lower fraction of 
warm gas at 2: = 0. The trend of decreasing WHIM 
fraction with resolution seen in the adiabatic runs has 
reversed. For z > 1, the fraction of hot gas increases 
with resolution. However, by z = 0, the difference in 
hot-gas fraction is negligible. As in the adiabatic simu- 
lations, the mass in the condensed phase increases with 
redshift, although the precise evolution with redshift is 
rather different. In general, resolution effects produce 
relatively minor changes in the evolution of the differ- 
ent phases. The most severe lack of convergence ex- 
ists for the warm and WHIM phases in the local feed- 



back simulations. The percent differences, defined as 
100 X |xi - X2I / {xi + X2)), in the warm and WHIM 
phases at z = between the highest and second-highest 
resolution simulations with 25 Mpc/h box size and local 
feedback are 12% and 17%, respectively. For the cor- 
responding distributed-feedback simulations, the percent 
differences for the warm and WHIM phases are 9.5% and 
4.4%. Thus, it is difficult to argue that the 25 Mpc/h 
simulations are completely converged. 

The larger box simulations (Figure [s]) show a greater 
degree of convergence. The percent differences for the 
50 Mpc/h simulations analagous to those given above 
are 6.9% (warm) and 4.0% (WHIM) with local feedback 
and 5.0% (warm) and 0.2% (WHIM) with distributed 
feedback. Since the star formation rate is the aspect of 
these simulations that is furthest from convergence, this 
may suggest that creation of WHIM gas via gravitational 
shock heating dominates over stellar feedback for larger 
box sizes. Unfortunately, due to computational cost, we 
were unable to investigate this by running any simula- 
tions with / > 50 Mpc/h at the same resolution as our 
most refined 25 Mpc/h simulation. However, as we con- 
tinue to optimize our numerical methods and the avail- 
able computational resources increase, we will perform 
additional simulations with larger box sizes and higher 
dynamic range. 

3.3.2. Redshift Evolution and Feedback 

All of our simulations display a behavior, seemingly 
unique to this work, with respect to the evolution of the 
WHIM gas fraction. The amount of WHIM peaks at 
z ^ 0.5 and is in decline when the simulation reaches the 
current era {z = 0). This coincides with the fraction of 
warm gas levelling out, as opposed to continuing to de- 
crease. Although there is no reason to believe that z = 
is a "privileged epoch," our result is in contrast to most 
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Figure 6. Evolution of baryon fraction in each phase for the three runs with 25 Mpc/h box size and 768^ grid cells. Shown are two 
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distributed feedback and y = 0.025 (dot-dashed). 



1.0 
0.8 
0.6 
0.4 
0.2 
0.0 



Warm 







1 2 3 4 
z 



0.5 
0.4 
0.3 
0.2 
0.1 
0.0 



WHIM 





\ 1 \ 
















\\ 




V 




1 







1 2 3 4 
z 



0.030 
0.025 
0.020 
0.015 
0.010 
0.005 
0.000 



Hot 



Condensed 



50_ 


.768_ 


_1 


- - 50_ 


.768_ 


.lb 


- ■ - 50_ 


_768_ 


.2 


50_ 


_768_ 


.2 b 
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simulations with local feedback, one with y — 0.025 (solid) and one with metal yield y = 0.005 (dashed), plus two simulations with 
distributed feedback, one with y = 0.025 (dot-dashed) and one with y = 0.005 (dot). 



other numerical studies, which find the WHIM fraction 




the WHlM traction to flatten out at low redshift rather 
th an cont inuing to increase. The "constant-wind" model 
of [Dave et al.j (|2010|) shows a decline in WHIM fraction 
at late time, although it is the only one of their sim- 
ulations to display thi s behavior. The "BH" model of 



Tornatore et al.| ([2010| also shows a decreasing WHIM 
fraction at low redshitt. However, that simulation, which 
was also the only one in that work to feature this result, 
was designed specifically to have feedback from super 
massive black holes, which we do not include. Hence, 
it is somewhat unclear why our results agree. Although 
the adiabatic simulations show this behavior as well, it 
is more pronounced in the simulations with cooling and 
feedback. 

In order to better understand this phenomenon, we at- 
tempt to measure the "flux" of gas from one phase to 
another in a rather crude manner. We consider flux be- 
tween phase A to phase B to exist when a single grid 
cell changes from phase A to phase B when measured 
in two temporally nearby data outputs of a simulation. 
The value of the flux from one phase to another is sim- 
ply taken to be the mass in the later data output of the 
cell that has changed phase divided by the time interval 
between data outputs. Unfortunately, we are unable to 
account for the change in amount of baryons in a given 
phase resulting from a change in density in cells that do 
not move from one phase to another. While this method 
is clearly flawed, to the extent that the exact numerical 
values should not be trusted, it does add useful qualita- 
tive insight. 

To disentangle the contributions of cooling and feed- 



back from cosmological evolution, we first perform this 
analysis on run 50_1024_0, the largest of our adiabatic 
simuations. For this exercise only, we divide the WHIM 
phase into two phases, the warm- WHIM (10^ K < T < 
10^ K) and the hot- WHIM (10^ K < T < 10^ K). As the 
fraction of gas in the warm phase decreases {z > 1), we 
measure a net flux from the warm phase to the warm- 
WHIM phase. At the same time, we also see a net flux 
from the warm- WHIM to the hot- WHIM. At low redshift 
{z < 1), where the evolution of the warm fraction ap- 
pears to flatten out, the flux reverses with warm- WHIM 
gas now moving back into the warm phase. Likewise, 
hot- WHIM gas now begins flowing back into the warm- 
WHIM phase. In simulations without radiative cooling, 
the only means of cooling is adiabatic expansion. This 
seems to suggest that, irrespective of cooling and feed- 
back, the universe undergoes an epoch of heating fueled 
by shocks from structure formation followed by an epoch 
of cooling due simply to the acceleration of the Hubble 
flow. For the cosmological parameters used in this work, 

1/2 

the expansion factor {E{z) = [^^^(1 + z)^ + I^a] ^ ) is 
dominated by for z < 0.4, which is consistent with 
the epoch at which the WHIM fraction begins to decrease 
in our simulations. When analyzing run 50_1024_2 in the 
same manner, we flnd a similar picture. In Figure [Sj 
we plot the net phase flux for the warm, warm- WHIM, 
hot- WHIM, and hot phases for run 50_1024_2. From this 
flgure, the epochs of heating and cooling are quite visi- 
ble, with the reversal of net input to net output occuring 
earlier in the higher temperature phases. In Figure [Ol 
we plot the contribution of each phase to the input and 
output of the warm and warm- WHIM phases. It is clear 
from Figure |9] that the warm- WHIM is supplied with ma- 
terial from cooler phases at early times duing the epoch 
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Figure 8. The net flux of baryons into and out of each phase in run 50_1024_2. Top-left: warm (T < 10^ K), top-right: warm- WHIM 
(10^ K < T < 10^ K), bottom-left: hot- WHIM (10^ K < T < 10^ K), bottom-right: hot (T > lO'^ K). 



of heating, and from hotter phases at late times during 
the epoch of coohng. Likewise, material leaves the warm- 
WHIM phase moving mostly into hotter phases earlier 
on, and then cooler phases later. 

The overall contribution of radiative cooling and feed- 
back to the evolution of the baryon phases is less clear, 
as they add both heating and cooling. With cooling and 
feedback included, the peak in the WHIM fraction occurs 
slightly earlier. In Figures [6] and [7| we focus on the influ- 
ence of feedback on the evolution of the baryon phases. 
For both box sizes, the simulation with distributed feed- 
back has roughly 10% more WHIM than the simulation 
with local feedback. The additional WHIM in the dis- 
tributed feedback simulations comes at the expense of 
condensed gas. This is another example of the ability 
of the distributed method to overcome the over-cooling 
problem associated with local feedback. Interestingly, 
the runs with lower metal yield, denoted with the letter 
"6^^ in their labels, reach their peak WHIM fraction at 
systematically higher redshifts than their higher metal- 
yield counterparts. The WHIM fraction peaks at z ~ 0.5 
in the higher-yield runs, but at z ~ 0.9 in the lower-yield 
runs, an offset of ~ 2.3 Gyr. At higher redshift, the low- 
yield runs have slightly more WHIM gas. However, the 
WHIM begins to decrease at earlier times, and the result- 
ing WHIM fraction at z = is lower. The low-yield runs 
also have lower condensed fractions, except in the case of 
run 50_768_2b. In every case, these offsets are balanced 
by a higher fraction of warm gas. In comparison to their 
high-yield analogs, the low-yield runs reach their peak 



SFR slightly earlier (2: ~ 2 as opposed to 2: ~ 1.5) and 
have significantly lower SFR at low redshift (roughly 0.5 
dex compared to their higher-yield counterparts). Both 
of these can be attributed to the lower radiative cool- 
ing efficiency resulting from the decreased metal yield. 
These clues highlight the importance of accurate treat- 
ments of cooling and feedback in the evolution of the 
WHIM. The use of distributed feedback increases the 
WHIM fraction by transporting stellar feedback into the 
IGM with greater efficiency. Lowering the metal yield 
produces a similar effect, but it also suppresses the SFR 
more effectively. Interestingly, the single simulation of 
Tornatore et al. (2010) to show a decline in WHIM frac- 
tion at low redshiff^rso happens to be the run in which 
the star formation rate is reduced by the greatest amount 
during that epoch. While radiative cooling may be re- 
sponsible for some additional amount of gas that is just 
above 10^ K, cooling into the warm phase, the majority of 
baryons residing in the WHIM are at low enough density 
that their cooling time exceeds the Hubble time. Thus, it 
seems that the role of radiative cooling is to enhance the 
ability of dense gas to cool and continue to form stars, 
which then inject thermal energy into the IGM, heating 
the gas into the WHIM phase. 

3.4. VI Absorption 

Ultraviolet absorption lines of O VI (1032, 1038 A) are 
believed to be excellent probes of the cooler portion of 
the WHIM, as the peak fraction of collisionally ionized 
O VI occurs at T ~ 3 X 10^ K. However,,Oppenheimer fc 
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Figure 9. Inputs (left panels) and outputs (right panels) for the warm phase (top panels) and warm- WHIM phase (bottom panels) in 
run 5U-1U24_2. 



[Dave ( |20Q9[ ) have claimed that the majority of O VI ab- 
sorption comes from warm, photoionized gas (T ~ 15,000 
K), although their results are likely to be affected by the 
significant over-cooling within their simulations. They 
suggest that O VI may not be a reliable WHIM tracer. 
Nonetheless, there exists a number of large s urveys of 



O VI absorption lines at low redsh ift (Danforth fc Shull 
2005; Danforth et al. 2006; Danfor th ^hull|:^00^| |'17ipp 



et al.,|2QQ8^ |Thom 
lent tests or numerical simulations, 
make such comparisons. 



Chen||2QQ8ajb ) that provide excel- 



In this section, we 



3.4.1. Creating Synthetic QSO Sight Lines 

Our simulations provide a useful database that can be 
compared to observations of O VI and other ions along 
sight line s to ac tive galactic nuclei (AGN). We use the 
Y'j^( Turk|2QQ8j [Turk et al. 2Q11J analysis toolkit to con- 
struct a set of 50U random AGN sight lines through each 
of our simulations. In ord er to compare with the O VI 
absorbers in the survey of Danforth & Shull (2008), we 
design each sight line to extend from z = to z = 0.4, 
the depth of the survey. Since the Sz corresponding 
to the length of our simulation boxes is much smaller 
than the desired redshift depth, our sight-line generation 
technique involves stacking multiple rays cast through 

^ http://yt.enzotools.org 



simulation boxes from different epochs. Before running 
the simulations, we compute the exact number of data 
outputs and their location in redshift space required to 
traverse the comoving distance from z = 0.4 to z = 0, 
using ray segments no longer than the box length. For 
example, at z = 0.4, a comoving distance of 50 Mpc/h 
corresponds to dz ~ 0.02. Thus, in this case, the next 
data output in the series is chosen to be at z = 0.38. We 
continue this process until reaching z = 0. For the 25 
(50) Mpc/h boxes, 47 (23) data outputs are required to 
span the redshift range. In this method, each sight line is 
able to sample the simulation box at continuously earlier 
epochs as it extends further away. Each individual sub- 
segment of a given sight line has a random (x, z) start- 
ing position and (^, (j)) trajectory so as to minimize the 
probability of sampling the same structures more than 
once per sight line. 

For each grid cell that is intersected by a ray, we 
record all relevant information of the grid cell, such as 
the baryon density; temperature; metallicity; redshift; 
path length, dl^ of the ray through the cell; and the cor- 
responding 5z. Each ray length element, or lixel, is as- 
signed a negatively incrementing redshift such that the 
first lixel has the same redshift as the associated data 
output and the last lixel has a redshift exactly 5z higher 
than the redshift of the next data output to be used. 
Note that since the ray segments are randomly oriented 
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Figure 10. Left: O VI ionization fraction, /ovi? in collisional ionization equilibrium (CIE) using data of [Sutherland Dopita| 
Right : /ovi calculated with Cloudy including collisional ionization and photoionization from a UV background aX, z = U Jjlaarcl 
|20QT] >. 



3 



and not necessarily parallel to any axis of the grid co- 
ordinate system, dl is not constant over an entire ray. 
The column density of the i'th ionic species of a given el- 
ement, X, is then given by Nx^ = x (^x/^h) x/xi x<^^ 
where (nx/nn) is the elemental abundance of X and 
is the ion fraction, or {ny^Jnx)- A coherent structure 
that may produce a single absorption line in reality is al- 
most certainly resolved over multiple lixels. In order to 
minimize cutting structures into multiple absorbers, we 
smooth each ray to a constant spectral resolution (A/(5A) 
of 5000, which is lower than the effective spectral reso- 
lution of a single grid cell. Quantities such as the tem- 
perature and metallicity of a smoothed lixel are taken as 
the mass- weighted averages of all contributing lixels. We 
do not find the resulting absorbers to be significantly af- 
fected by changes of up to an order of magnitude in either 
direction to the chosen spectral resolution. 

Since we do not track individual heavy element species, 
we assume the metals are present in solar abundances 
patterns such that (nx/^n) = Z x (nx/^n)©- For the 
oxygen abundance, we use the value of (no/nn)© = 
4.9 X 10~^, as recommended by Asplund et al.| (|2005|), 
and we employ two different methods to calculate the 
ionization fraction of O VI. We first consider the extreme 
case of no ionizing radiation, where O VI is in collisional 
ionization equilibrium (CIE) and /ovi is only a func- 
ti on of temperature. For this, we use the calculations 
of Sutherland & Dopita (|1993). We note the existence 
of new er calculations performed by Gnat fc Sternberg 
(2007), based on updated ionization and dielectronic re- 
combination rates. We intend to use these newer rates 
in our future work. Fortunately, the values for O VI dif- 
fer only marginally between the two studies. We plot 
/ovi(T) for CIE in the left panel of Figure [To] In the 
second method, we assume t he sam e UV metagalactic 



background (Haardt fc Madau|2001 ) used in the simula- 
tions. We do not apply any additional renormalization 
to the intensity of the spectra. At 2; = 0, the mean in- 
tensity, Jy^2ii of the background in units of 10~^^ erg 
cm-2 s-^ ilz-^ is 3.8 x 10-^ at 13.6 eV and 2.6 x lO""^ 
at 138 eV (the ionization energy of O VI). At z = 0.4, 
J^,2i(13.6 eV) = 0.12 and J^,2i(138 eV) = 7.3 x 10"^. 
When photoionization is included, the equilibrium value 
of /ovi becomes a function of density as well as tempera- 
ture. In the right panel of Figure 10, we show /ovi(p, T) 



calculated with Cloudy for the UV background at z = 0. 
For proper number densities uh ^ 10~^ cm~^, the equi- 
librium value of /ovi is determined solely by electron 
collisions. At lower densities, photoionization begins to 
dominate, and the metals are typically in a higher ion- 
ization state at a given temperature than in CIE. Thus, 
the temperature of maximum O VI fraction decreases 
with decreasing density. There is significant evol ution in 
the int ensity of the UV background of Haardt fc Madau 
( |2001[ ) between z = 1 and z = 0. We account for this 
when constructing the synthetic AGN sight lines by cal- 
culating /ovi values for each pixel intersected by the ray 
segment, using interpolation tables that vary in density, 
temperature, and redshift. Hereafter, we will refer to this 
method for calculating /ovi as the collisional + photoion- 
ization (C+P) model. 

3.4.2. Redshift Distribution of O VI Absorbers 

Numerous studies have attempted to recreate the ob- 
served number density per unit redshift of O VI absorbers 
as a function of column density to gauge the accuracy of 
their simu lations jFang & Bryan 2001 ; Cen & Fang 20Qi6l 
Oppenheimer & Dave 2009; Te pper-Garcia et al. ||20 10[ ) . 
In Figure [11] we plot the distribution of U VI absorbers 
per unit redshift per log A^ovi {d'^M / dzd(\ogNowi)) for 
the simulations with 25 Mpc//i box size and metal yield, 
y = 0.025. We separate runs by feedback method and 
plot O VI distributions calculated from both the CIE 
and C+P models. The number of O VI absorbers in- 
creases systematically with increasing resolution. This 
is not surprising, as the star formation rate, and hence 
the production of metals, also incr eases with resolution. 
Howev er, it should be noted that "Tep per-Garcia et al. 
(2010) show the opposite trend (top-left panel of Figure 
B2 in Appendix B of that work), with higher resolution 
simulations producing fewer O VI absorbers. Just as 
in the evolution of the baryon phases, the simulations 
with distributed feedback are more converged than those 
with local feedback. The number of O VI absorbers at 
very high column density (A/'ovi ^ 10^^ cm~^) is consis- 
tent at all resolutions. This is not surprising, since these 
absorbers are probably associated with large, collapsed 
structures that even the coarsest simulations are capa- 
ble of resolving. The disagreement between simulations 
with different resolution is greater for lower column den- 
sity. In all cases, the lowest resolution simulations are 
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Figure 11. Convergence of d'^Af/dzd{\og Noyi), the number density of O VI absorbers per unit redshift and per A(\ogNoyi), for the 
simulations with 25 Mpc//i box size and metal yield, y = 0.025. Line colors denote the resolution of the simulations as in Figure |2] The 
top two panels show the simulations with local feedback, while the bottom panels show those with distributed feedback. The left panels 
represent O VI fractions calculated with the CIE model, while the right panels represent O VI fractions calculated with the C+P model. 



grossly under-resolved with respect to the production 
of low column density absorbers, where we do not see 
full convergence. For the three highest resolution sim- 
ulations, the increase in the number of O VI absorbers 
is roughly constant with increasing resolution. It is not 
completely clear what resolution is required to reach con- 
vergence. However, as we discuss below, we are able to 
achieve an excellent match with the observed distribution 
of O VI absorbers with our highest resolution distributed- 
feedback simulation using the C+P method, which also 
happens to show the greatest convergence. 

In Figures [12] and lis) we plot the cumulative number of 
O VI absorbers, dJv{>Noyi)/dz^ given by the integral of 
the bivariate distribution {d'^M / dz dNo^i) over column 
density A/'ovb for our highest-resolution simulations with 
both 25 and 50 Mpc//i boxes. For all simulations shown, 
we plot the dM / dz resulting from /ovi calculations, us- 
ing both the CIE model (left panels) and the C+P model 
(right panels). We compare our values with the observa- 
tions of panforth fc Shull| (|2QQ8). The C+P model pro- 
vides a much better tit to tne"observed dj\f / dz for A^ovi 
< 10^^ cm~^. In contrast, the CIE model produces an ex- 
ces s of lo w column density O VI absorbers. As we discuss 
in §3.4.4[ these absorbers exist primarily in regions of low 
iryon d* 



local feedbac k. The semi-analytical models of Furlanetto 
et al. ( 2QQ5[ ) show steep drop-off in the number of O VI 
absorbers at A^ovi ^ few x 10^^ cm~^, which is approxi- 
mately where we observe dM{>NoYi )/dz t o drop below 
unity nearly universally. We show in §3.4.5| that the local 
feedback simulations have significantly more metal in re- 
gions of high overdensity, which is likely responsible for 
the increase in high column density absorbers. For both 
feedback methods, the CIE model is in closer agreement 
with observations for the high column density absorbers. 
An explanation for this is not immediately clear and may 
require the use of full radiation transport hydrodynamic 
simulations instead of the spatially uniform background 
employed in this work. 

In general, the simulations with the standard metal 
yield, y = 0.025, display broadly similar results, indepen- 
dent of feedback method. The runs using the distributed 
method show a higher dM /dz for A^ovi ^ 10^^ cm~^. 
This difference is more significant when using the CIE 
model. The increase in low column density absorbers 
from local to distribu ted feedback niodels was also seen 
in the simulations of Cen & Fang (2006). When pho- 
toionization is included in the /ovi calculation, the runs 
with local feedback produce more O VI absorbers for 



baryon density and are more infiuenced by ionizing radi- A^ovi ^ 10 cm ^ than the runs with distributed feed- 



ation. For higher column densities (A/'ovi ^ 10^ cm~ ), 
the C+P model creates slightly too many absorbers. 
This is within the error bars for the simulations with dis- 
tributed feedback, but somewhat outside for those with 



back. With the CIE model, there is essentially no dif- 
ference in the number of high column density absorbers 
between the two feedback methods. The dM /dz values 
for the 50 Mpc//i boxes show behavior similar to those 
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Figure 12. Cumulative number density of O VI absorbers per unit redshift as a function of column density for 500 random sight lines 
extending from z = to z = 0.4 for runs 25_768_1, 25_768_lb, and 25_768_2, compared with the observational data of Danforth &; Shull] 
(12008). Left: /ovi calculated assuming collisional ionization equilibrium (CIE). Right : fn\n calculated includ ing both collisional ionization 
and photoionization (C+P) from a redshift-dependent UV metagalactic background ( jHaardt Madau|200l| >. 




log Novi [cm ^] 

Figure 13. Same as Figure except for runs 50_1024_1 and 50_1024_2. 

of the 25 Mpc/h boxes, although they are shghtly lower 
in general. This is most likely due to the fact that the 50 
Mpc/h simulations are slightly less converged than the 
25 Mpc/h simulations. However, the dJV/dz values do 
not appear to show any specific dependence on box size. 

As expected, lowering the metal yield decreases the 
number of absorbers at all column densities. However, 
the decrease in absorber numbers is not constant over 
the range of column density. For A/'ovi > 10^^ cm~^, 
both the CIE and C+P models produce roughly 3 times 
fewer O VI absorbers in the low-yield run. In the CIE 
model, this ratio increases to a peak value of '^12 at 
^ovi ^ 10^^-^ cm~^, where it remains constant up to the 
maximum column density in the sample. For the C+P 
model, this ratio continues to increase over the entire 
range of column densities to a value of ~40 at A/max- 

3.4.3. Qovi 

The ratio of the total O VI density to the closure den- 
sity of the universe, or ftgyj, has been calculated for 
high CCarsweh et al.|l2002| jSimc oe et a'L||2QQ 4l |Ber^eron| 



log Nqvi [cm ^] 



Table 2 

^7ovI(xlo-7) 



Run 


Simulation Box 


Absorption Lines 




at z 


= 


log(A^/cm 


-2) > 12 (13) 




CIE 


C+P 


CIE 


C+P 


25_768_1 


17 


11 


7.8 (4.3) 


4.8 (4.5) 


25_768_2 


20 


3.5 


15 (9.4) 


3.0 (2.5) 


50_1024_1 


7.3 


3.4 


5.1 (2.9) 


3.0 (2.8) 


50_1024_2 


10 


2.5 


7.2 (4.5) 


2.3 (2.1) 



Note. — The contribution of O VI to the closure 
density of the universe for the four highest resolution 
simulations with normal metal yield. In columns 2 and 
3, r^ovi is calculated by summing the total mass in O VI 
in the simulation box at z = using the CIE (column 
2) and C+P (column 3) models. In columns 4 and 5, 
^OVI is calculated from our collection o f synth etic O VI 
absorbers via the method described in §3 4.3| 



& Herbert-Fort 2005) and low ( |'l'ripp et al. 



000; Dan- 

Ibrth & Shull 2008) redshift from the observed density of 



O VI absorption systems. Since all of these surveys have 
an inherent minimum column density to which they are 
sensitive, an important question is how much O VI goes 



undetected. Because the column densities of our syn- 
thetic O VI absorbers are calculated directly from the 
simulation data and not derived from analysis of spectra, 
our synthetic survey has essentially infinite sensitivity, 
limited only by the spatial resolution of the simulations. 

We calculate ftoyi from our synthetic O VI survey by 
integrating over the number density of absorbers per cos- 



Nature of the WHIM 



15 



CIE C+P 



25_768_1 




Figure 14. l^ovi as a function of minimum O VI column density for 500 random sight lines extending from 2; = to 2; = 0.4 for runs 
25_768_1, 25_768_lb, and 25_768_2, compared with the observational data of Danforth & Shull (2008). Left: /ovi calculated assuming 
collisional ionization equilibrium (CIE). Right: fovi calculated includin g both collisional ionization and photoionization (C+P) from a 
redshift-dependent UV metagalactic background ( [Haardt Madau|2001| ). 
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Figure 15. Same as Figure [M] except for runs 50_1024_1 and 50_1024_2. 



mological path length 
^0^0 VI 
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This is expressed as 
dXdlogN 



Jo JNmi 



NdlogNdX, 
(4) 

where Hq is the Hubble constant in units of s ^, movi 
is the O VI mass, c is the speed of light, pcr,o is the 
critical density of the universe at z = 0, and dX is th e 
''absorp t ion pa th length function" , defined by [Bahcall 
Peebles! (|1969D as 



-1/2 



dz. 



(5) 



The absorption path length function corrects the number 
density of absorbers per unit redshift for the expansion 
of the universe. For the redshift interval from z = to 
0.4, X 0.53. Since all of our synthetic sight lines probe 
this redshift epoch, the total path length used in Equa- 
tion H] is simply 0.53 multiplied by the total number of 



sight lines. In Figures [M] and [15) we plot l^ovi as a func- 
tion of the minimum O^VI column density. The values 
of r^ovi are reflective of the values of dM /dz and, there- 
fore, the results are very similar. The best matches to 



the observational data of Danforth & Shull ( 2008 ) come 
from the simulations with distributed feedback and using 
the C+P model for calculating /ovi- Due to the shal- 



low slope of dM /dz for low column density absorbers in 
the C+P model, the additional O VI mass below current 
sensitivity limits of A^qvt ^ 10^^ cm ~^ appears to be neg- 
ligible. However, [Danforth fc ShuH] ( [2008| ) find the differ- 
ential column density distribution (WTT/ dzdilogNoYi)) 
to be proportional to ~ A^~^, implying equal mass at all 
column densities. Determination of the slope of the O VI 
distribution at low column density is a primary science 
goal for the newly installed Cosmic Origins Spectrogaph 
on HST. 

We have the luxury of also calculating Q.0VI via direct 
summation of all the grid cells in the simulation box. 
This allows us to determine whether calculation of O^ovi 
in the above manner produces a biased result. In Table 
[2) we list the values of Q^ovi calculated through direc- 
tion summation using both the CIE and C+P models 
for each of our four highest resolution simulations. We 
also give the values calculated from our synthetic ab- 
sorption line samples for minimum column densities of 
10^^ cm~^ and 10^^ cm~^. The two methods are gener- 
ally in good agreement, with direct summation yielding 
a value roughly 10-40% higher, with the exception of run 
25_768_1 in which both the CIE and C+P summed values 
are approximately double their counterparts. 
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Figure 16. Mean H number density, metallicity, temperature and 
O VI fraction at z = 0, for two methods of feedback distribution. 
Points show 1 a devi ations, with column densities of O VI absorbers 
shown in Figure [12] (25 Mpc//i boxes with metal yield y = 0.025). 
Black lines show statistics for O VI fractions using CIE model and 
grey lines represent O VI fractions calculated with C+P model. 

3.4.4. The Physical Conditions of O VI Absorbers 



In Figures 16 and [TtI we examine the physical condi- 
tions of the 0> VI absorbers plotted in Figures [12] and 
13 The choice of /ovi ionization model has consider- 
able influence on the average physical conditions of the 
O VI absorbers as a function of their column density. 
For both the CIE and C+P models, the average gas den- 
sity per absorber rises with column density. For the CIE 
model, we find uh oc A/'qyj, and for the C+P model, we 
find uh cx A/'qyj. The slope of this increase is signif- 
icantly steeper with the CIE model. At the minimum 
column density measured, the average gas density as- 
sociated with O VI absorbers from the CIE model is 
approximately 1.5 orders of magnitude lower than with 
the C+P model. For A/'ovi > 10^^ cm~^, the associ- 
ated gas density is roughly similar for both models, with 
UH ^4..5x 10-^ cm-3 or 5h ^ 240. With WMAP-7 pa- 
rameters, Qbh^ = 0.02260 ± 0.00053, the baryon density 
and corresponding hydrogen number density uh can 
be written 



Pb-- 



■- (4.26 X 10"^^ g cm-^)(l + zf 5h 



Uh = (1.90 X lO"'^ cm-^)(l + zf Sh , 



(6) 



where Sh is the overdensity factor. 

The spatial difference between O VI in the CIE model 
and the C+P model is made clear in Figure [18] where we 
show projections through the full simulation box at z = 
of mass- weighted mean baryon density and temperature, 
along with O VI column density, using both /ovi models 
for run 50_1024_2. With the CIE model, O VI exists spa- 
tially much further from collapsed structures, where the 
gas temperature is in the optimal range for O VI. How- 
ever, with the C+P model, gas that would have a high 
/ovi in CIE is instead photoionized to higher ionization 
states. In this case, O VI is confined primarily to IGM 
filaments. 
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Figure 17. Same as F igu re [Tg] but for the runs with the 50 Mpc/h 
boxes shown in Figure [Ts] 



With the CIE model, the metallicity of O VI absorbers 
is relatively independent of column density and consis- 
tent with Z = Zq up to A^ovi ^ 10^^ cm-^. With 
the C+P model the metallicity rises with column den- 
sity from - 10"^ Zq at TVqvi = 10^^ cm-^ to Zq at 
10^^ cm~^. The C+P model results are similar to the 
findings of [Cen fc Fang (2006), who see a rise from just 
over 0.1 Zq to slightly under Zq. The density and metal- 
licity values from th e C+P models are broad l y con sistent 
with the results of Oppenheimer & Dave ( 20091 ), who 
employ a similar method lor calculating /ovi- Ihe av- 
erage absorber temperature remains fairly constant as a 
function of column for both /ovi niodels. For the CIE 
model, the average temperature is ~ 3 x 10^ K, with a 
small variance. This is not surprising, as the peak value 
of /ovi in CIE is sharply peaked around this temper- 
ature. For the C+P model, the average absorber tem- 
perature is ~ 10^ K. Unlike the density and metallicity 
values, this is in consider ably less agreement with |0p-| 
penheimer & Dave (2009), who find the mean absorber 
temperature to be ^ 15,000 K. This is most likely due 
to the higher cooling rates in that work resulting from 
not accounting for the radiation background when com- 
puting the cooling from metals. The mean value of /ovi 
is roughly 0.1 with the CIE model over the entire range 
of column densities. The values of /ovi from the C+P 
model are typically lower than the CIE values by ~ 0.5 
dex for A/'ovi ^ 10^^ cm~^. The highest column density 

1015.5 



absorbers (A'ovi 



> 



cm ^) tend to exist as a re- 



sult of particularly optimal conditions for O VI: high gas 
densities in excess of 10~^ cm~^, super-solar metallici- 
ties, and values of /ovi near the maximum allowable in 
either CIE or C+P models. 

For the highest resolution 25 Mpc/h simulations with 
normal metal yield, the run with local feedback forms 
approximately 66% more stars, and has essentially the 
same proportion of metals, compared to the run with 
distributed feedback, yet it produces 38% fewer O VI 
absorbers with A^ovi ^ 10^^ cm~^. How is this pos- 
sible? To answer this question, we examine the distri- 
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Figure 18. Projections through the full simulation box of run 50_1024_2 at 2; = of mass- weighted mean baryon overdensity (top- left), 
mass- weighted mean temperature (top-right), and O VI number density, with /ovi calculated assuming CIE (bottom- left) and including 
photoionization (C+P, bottom-right). Projections are constructed by summing the values in all grid cells along the line of sight. For 
weighted projections, each pixel in the image is calculated as ^rriiWi/ ^Wi^ where rrii is the projected quantity and Wi is the weighting 

i i 

quantity. For more information, see |Turk et al.| ( [2011| ). 

bution of gas-phase baryon (Figure p^O]) and metal mass 
(Figure 20) in bins of baryon overdensity and temper- 
ature for runs 25_768_1 and 25_768_2. In runs 25_768_1 
and 25-768-2, ~ 75% and 85% of the total baryon mass is 
in the gas-phase, as opposed to in stars. The distribution 
of baryons in the two simulations is qualitatively similar, 
but there are subtle differences. The simulation with lo- 
cal feedback has a significant population of warm, high 



density baryons (T < 10^ K, pb/pb ^ 10 ) 5 ^^t present 
in the simulation with distributed feedback. Instead, 
the distributed-feedback simulation shows an increased 
amount of hotter, slightly less dense baryons (T ^ 10^ 
K, pb/pb ^ 10^~^), and more baryons in lower density 
WHIM as well. Recall from Figure [6] that the local feed- 
back run has a WHIM mass fraction of 34% at z = 0, 
compared to 44% in the distributed-feedback run. 

When considering the distribution of metals, the two 
simulations appear even more distinct. A significant por- 
tion of the gas-phase metals in the local feedback run ex- 
ist in the warm, high-density gas that is not present in the 



distributed-feedback run. These metals are newly cre- 
ated by stars, yet unable to escape their points of origin 
due to over-cooling. While the two runs produce roughly 
the same number of gas-phase metals, the run with lo- 
cal feedback produces significantly more metals in stars. 
Defining ^metals,* and ftmetais,gas to be the total amount 
of metals in stars and gas with respect to closure den- 
sity of the universe, we find ^metals,* = 1-17 x 10~^ and 
^metals, gas = 3.93 X 10~^ for the ruu with local feedback 

and ftmetals,* = 4.81x10""^ and ftmetals,gas = 3.89 xlO""^ 

for the run with distributed feedback. In contrast, the 
distributed-feedback simulation shows far more metals 
in the WHIM phase. To be exact, of the metals not in 
stars, 44% by mass are in the WHIM phase in the lo- 
cal simulation, whereas 73% are located in the WHIM in 
the distributed simulation. Differencing the two panels 
of Figure 20 reveals that the greatest increase in metal 
mass, from the local to the distributed simulations, lies 
on an adiabat extending from T ~ 10^ K and pb/pb ^ 10^ 
directly into the heart of the WHIM. These are metals 
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Figure 19. Fraction of total gas-phase baryon mass in two-dimensional bins of overdensity (pb/pb) and temperature for runs 25_768_1 
(left) and 25_768_2 (right) at 2; = 0. The concentration of baryon s in the thin red line extending from Pb/pb ^ 10~^ and T ~ 300 K up to 
Phi Pb ~ 10^ and T ~ 10^-^ K represents unenriched (see Figure [20] for an indication of metallicity) IGM lying on an adiabat. Radiative 
cooling becomes important at higher densities, as evidenced by the declining baryon temperature. The baryon concentration in the regime 
of 10^ ^ Pb/ Ph ^ 10^ and 10^ K < T < 10^ K depicts collapsing halos, galaxies, and star- forming regions. Virialization shocks and stellar 
feedback heat gas up to ~ 10^ K where it travels out into the IGM, cooling adiabatically. The adiabat extending from Pbl Pb ^ 10^ and 
T ~ 10^ down into the broad plume of baryons between 10^ K and 10^ K is the WHIM. The fraction of total baryon mass in gas (not 
stars) is 75% for run 25_768_1 and 85% for run 25_768_2. 
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Figure 20. Fraction of total gas-phase metal mass in two-dimensional bins of overdensity {pb/ Pb) and temperature for runs 25_768_1 (left) 
and 25_768_2 (right) at z = 0. The primary regimes of metal population are the galaxies and the WHIM. Although the total amount of 
gas-phase metals produced is roughly equivalent between the two runs, the amount of metals in stars is much higher for run 25_768_1. The 
fraction of the total metal mass in gas (not stars) is 25% for run 25_768_1 and 45% for run 25_768_2. 



outflowing from galaxies into the IGM, cooling adiabat- 
ically as they expand into under dense gas. In summary, 
even though the local feedback simulation produces more 
metals, most of them become locked up in their host 
galaxies, unable to escape into the IGM. In contrast, 
the distributed-feedback simulation creates overall fewer 
metals, but it transports them into the WHIM far more 
efficiently. 

3.4.5. Does O VI Trace the WHIM? 

In order to answer this question, we plot the distri- 
bution of O VI mass from the C+P model in bins of 
temperature and baryon overdensity for run 50_1024_2 



at z = in the left panel of Figure [2Tj From this, we can 
see that the impression given by Figure [16] that the mean 
temperature of O VI absorbers is 10 K is misleading. 
Instead, it appears that O VI has a bimodal distribution 
in temperature, with populations existing at ~ 10^*^ K 
(mostly photoionized) and at ~ 10^*^ K (primarily colli- 
sionally ionized). We do not see any significant fraction 
of O yi at T ~ 15, 000 K, as claimed by Oppenheime r fc| 
Dave ( 2009| ). For the simulation shown, the fractions of 



total O VI mass are: 55% in the WHIM, 37% in the warm 
phase, and 8% in the condensed phase. The fraction of 
O VI in the hot phase (T > 10^ K) is negligible. 
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Figure 21. Left: fraction of total O VI mass in two-dimensional bins of baryon overdensity (pb/Pb) cind temperature for run 50_1024_2 at 
z = 0. Right: baryon overdensity and temperature associated with each synthetic O VI absorber in run 50_1024_2 (right panel of Figure 
|13| . Colors correspond to the colu mn d ensity of the absorber. In each panel, solid lines delineate divisions between the warm, WHIM, and 
condensed phases, as described in §3.3| 

cell results in an unphysically high cooling rate. We find 
that using the distributed-feedback method significantly 
lowers the low-redshift SFR and provides a much bet- 
ter fit to observations. For both feedback models, the 
SFR appears converged with respect to spatial and mass 
resolution at low redshift. At high redshift, the simu- 
lations are far from converged, even though the SFRs 
of the highest-resolution runs are not in significant dis- 
agreement with observations, however unreliable these 
observations may currently be. Nevertheless, the final 
mass in stars is reasonably well converged, as the bulk of 
cosmic time occurs at low redshift. 

The choice of feedback method has a significant im- 
pact on the fraction of baryons in the WHIM. For our 
highest-resolution simulations, there is a 10% difference 
in WHIM fraction at z = 0, with ~ 35% of baryons in 
the WHIM in the local-feedback simulation and ~ 45% in 
the distributed-feedback simulation. However, these re- 
sults are not totally converged. There is a significant lack 
of convergence in the WHIM fraction for the 25 Mpc/h 
box simulations with local feedback. The distributed- 
feedback simulations appear much closer to convergence 
on the WHIM fraction. However, those runs still show 
a slight increase in condensed gas and a decrease in dis- 
tributed gas with increasing resolution. It is encourag- 
ing that the baryon fractions for the distributed-feedback 
simulations of both box sizes are in good agreement with 
each othe r. In slight contrast to other work on this sub- 



Given these statistics, the next question is: Are these 
phases proportionally represented by the VI absorbers? 
In the right panel of Figure [2l] we plot each one of the 
synthetic O VI absorbers created for run 50_1024_2, us- 
ing the C+P model, on the plane of baryon overdensity 
and temperature with which they are associated, coloring 
each absorber according to column density. With respect 
to the total column density, 59% of the synthetically ob- 
served O VI is in the WHIM, with 36% in the warm phase 
and 5% in the condensed phase. Given the somewhat ar- 
bitrary definition of the condensed phase, the calculated 
fraction of O VI in this phase is not likely to be a very 
robust quantity. If we consider only O VI absorbers with 
column densities greater than 10^^ cm~^, closer to the 
observable range of O VI with HST and FUSE, the frac- 
tion in the WHIM decreases to 57%, the warm fraction 
increases to 37%, and the condensed fraction increases 
to 6%. This indicates that there is little biasing of the 
observed O VI away from t he true O VI distribut i on. In 



(2010) 



contrast, the simulations of |Tepper-Garcia et al 
show a similar distribution of O VI mass (bottom panel 
of Figure 3 in that work), yet do not produce a significant 
number of O VI absorbers in the warm phase. 

4. DISCUSSION AND SUMMARY 

We have performed a series of cosmological simulations 
designed to study the evolution of the WHIM and its ob- 
servability in low-redshift O VI absorption lines. The 
primary goals of this work are to find an optimal param- 
eterization of the processes of star formation and feed- 
back and to understand the extent to which the results 
of these simulations are converged and can be trusted. 
We investigated two methods of depositing stellar feed- 
back onto the simulation grid: (1) injecting all of the 
gas, metal, and thermal energy created by a star particle 
into the single cell in which it is located (local feedback); 
(2) distributing this feedback evenly over the central cell 
and its 26 nearest neighbors (distributed feedback). At 
low redshift {z ^ 1), the local-feedback simulations show 
star formation rates (SFRs) much higher than observed. 
This appears to be the result of over-cooling that occurs 
when the injection of thermal energy into a single grid 



ject, e.g.. 



Dave et al. (2001), we 



|Cen fc Ostriker| ( 2006 ) ; 

do not see the traction oi WHIM gas continuing to rise 
at z = 0. Instead, all our simulations show the WHIM 
fraction to peak at z ~ 0.5, then decrease by a few pe r- 
cent by the present day. As noted above in §3.3.2[ a 
small number of simulations in other works have shown 
a decline in WHIM fraction at late time. However, none 
of those studies see this phenomenon to be as ubiqui- 
tous as we do. We find that the decline in WHIM is 
offset by a leveling off of the warm gas fraction that is 
decreasing until z ^ 0.5 and a continual increase in the 
condensed fraction. By measuring the fiux of baryons 
from one phase to another, we see that the universe un- 
dergoes an epoch of heating fueled by the formation of 
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structure, followed by an epoch of adiabatic cooling when 
the Hubble flow begins to accelerate due to the influence 
of Encouragingly, the redshift at which I^a becomes 
dominant {z ~ 0.4, with the cosmological parameters 
used here) is consistent with when we observe the WHIM 
fraction starting to decrease. Radiative cooling has two 
effects. The total amount of condensed gas is increased 
by the addition of radiative cooling, coming at the ex- 
pense of WHIM gas. However, radiative cooling, primar- 
ily from metals, also enables dense gas to recool after 
being heated by stellar feedback, allowing it to continue 
to form stars and inject heat into the IGM. 

We used our highest-resolution simulations to create 
populations of synthetic O VI absorbers. We compared 
these with the observed number density per unit redshift 
{dM/dz) of O VI absorbers from z = to z = 0.4 as 
measured by Danforth & Shuir(2008). We experimented 
with two methods for calculating the fraction of oxy- 
gen in O VI: the CIE model, which assumes collisional 
ionization equilibrium, and the C+P model, which adds 
liotoionization from the UV background of |Haard"t 



pi 

|m 



Madau| ( |2001^ ). Overall, the distributed-feedback simula- 
tions, m conjunction with the C+P model for comput- 
ing /ovb provide the best match to the observations. 
The C+P model provides an exceptional fit for low col- 
umn density absorbers with Nqyi < 10^^ cm~^, but it 
slightly over-predicts the number of higher column den- 
sity absorbers, although still within the error for the 
distributed-feedback simulations. The CIE model fits the 
dM/dz of the higher column density (A/'ovi ^ 10^^ cm~^) 
absorbers better than the C+P model. Previous stud- 
ies have also noted that low column density absorbers 
are better matched with photoionization models, while 
high column density absorbers ar e better matched with 
collis ional ion ization models (Cen et al. 2QQT| [Fang fc 



Bryan 2001; Chen et al.||2003| ). |l^'urlanetto et al.| ( |2005p 

used relatively simple arguments to predict the maxi- 
mum column density of characteristic O VI absorbers to 
be slightly over 10^^ cm~^. In nearly all of our models, 
we observe the dM{>NoYi)/dz of O VI to drop below 
unity at roughly similar column densities. Cen & Fang] 
|2006) noted, as do we, that there is some dependence 
ot the number of low column density absorbers on feed- 
back method, with distributed-feedback methods pro- 
ducing more than local- feedback methods. This differ- 
ence is much more significant when using the CIE model 
to calculate /ovi- As is shown in Figure [18) the CIE 
model produces large amounts of O VI at great distances 
from collapsed structures, precisely where distributed- 
feedback methods excel over local feedback. On the 
other hand, the C+P model limits O VI to regions much 
closer to galaxies where the choice of feedback model is 
less likely to make such a difference. Given that the 
C+P model appears to be a much better fit to the ob- 
served number density of low column density absorbers, 
the prospects for constraining the properties of galac- 
tic outflows with O VI may be limited. In ionization 
equilibrium with the C+P model, collisional ionization 



is dominant for n > 10 cm (Figure 10). However, 



the mean baryon density associated with OVI absorbers 
with the high column density O VI absorbers is between 
~ 4 X 10~^ cm~^ and ~ 4 x 10~^ cm~^, or overdensities 
ranging from 5h ^ 23-2,300 at z = 0. Future simula- 
tions that incorporate radiation transport will likely be 



required to fully understand why the influence of pho- 
toionization appears to be reduced in these regions. 

The average physical conditions sampled by the two 
O VI ionization models are quite different, especially for 
low column densities. The highest column density ab- 
sorbers, with A/'ovi ^ 10^^ cm~^, have similar properties 
in both models and typically result from the rare combi- 
nation of high physical density, high metallicity, and high 
/ovi- In general, O VI absorbers from the CIE model 
have higher associated baryon densities, lower metallic- 
ities, higher temperatures, and higher /ovi than those 
in the C+P model. In a later paper, we will investigate 
the implications of the covariance in these parameters 
for the (O Vl-based) IGM baryon census. While O VI 
absorbers from the CIE model probe temperatures in a 
narrow range around the peak in /ovi (~ 300, 000 K), the 
average O VI absorber temperature in the C+P model 
is much closer to 100, 000 K. This is noteworthy, as it 
differs significantly from t he recent simulations of |0p 



penheimer & Dave (2009) who find O VI absorbers to 
have average temperatures of ^ 15, 000 K. Although the 
mean absorber temperature is ~ 100, 000 K in the C+P 
model, a closer inspection reveals two distinct popula- 
tions of O VI at - 30, 000 K and - 300, 000 K. We find 
that 55% of ah O VI is in the WHIM, 37% is in warm gas, 
and 8% is in condensed gas. These proportions, which 
come from summing the total O VI mass in the simu- 
lation box at z = 0, are reasonably well represented by 
the sample of synthetic O VI absorbers, with only a few 
percent bias toward the WHIM. 

We acknowledge a number of caveats to our current 
work. Radiation transport and ionizing sources are not 
explicitly present in our simulations, nor do we include 
any ionizing sources or spatially varying radiation back- 
ground. The radiation background is treated by alter- 
ing the H/He photoionization rates in a time varying, 
yet spatially uniform manner. The addition of radiation 
transport and ionizing sources may change the results 
significantly, given that the O VI absorbers are observed 
to be clustered ar ound filaments (Stocke et al. 2006) 
Wakker fc Savage] [200 9 ) and thus statistically nearer 
to sources of ionizing radiation. We also do not track 
the no n-equilibr i um ev oluti on of /ov i, as has been done 
by ,Cen fc Fang ' (2006) and [Yoshikawa fc Sasakij ( |2006| ). 
While Cen fc Fang (2066|) showed that non-equilibrium 
calculations of /ovi can alter the dM/dz distribution in 
O VI column density, we were still able to fit the obser- 
vations remarkably well with the C+P model. In future 
work, we will include a non-equilibrium calculation of 
the evolution of ionization fractions for oxygen and other 
heavy elements. 

We have shown that the number of O VI absorbers de- 
pends sensitively on the metal yield from stars, but we 
have not included feedback processes that rel ease metals 
in varying abundanc e patterns. For instance, Dalla Vec 



chia fc Schaye ( 2008 ) consider the feedback from Type la 
and 'iype II supernovae separately. This also changes the 
rate of injection of thermal energy. The parameters we 
have chosen for our star formation and feedback methods 
reflect our assumption that stars form with a Salpeter 
IMF at all ti mes, which is not necessarily true. 



et al. (2008) have pointed out inconsistencies between 
the observea star formation rate and the observed stellar 
density; the observed SFR implies a much higher stellar 



Nature of the WHIM 



21 



density than observed. Since we are able to match the 
observed evolution of SF R, our simula t ions a re also pro- 
ducing too many stars. [Wilkins et al. (|2008') claim that 
this problem could be solved if the ste 



lar IMF was more 
to p-heavy at high redshift ^ 3), an idea also supported 
by Tumlinson (2007). Other recent works have cited ob- 
servational discrepancies that may be potentially solved 
with an evolving IMF with characteristic mass i ncreas- 
ing with reds hift, e.g.,| Davel ( |2QQ8| ; |van Dokkum| ([2008 ) . 
[Smith et al.| ( [200 9) have shown that the high-redshift 
TJMB can suppress fragmentation in star forming clouds, 
leading to an IMF that evolves with cosmic time and cre- 
ates additional massive stars in the past. A varying IMF 
would have a significant effect on the creation and disper- 
sal of metals in the early universe, and may be sufficiently 
important to consider in studies of the enrichment of the 
IGM. 

In subsequent papers, we will address the observability 
of the WHIM in O VII and O VIII and other important 
elements, such as C, N, and Ne. We will also investigate 
other radiation backgrounds and their influence on pho- 
toionized absorbers. The slope of dM /dz at low column 
density is sensitive to the shape and intensity of the radi- 
ation background, offering a means of constraining this 
background with simulations. In the very near future, 
it will be feasible to include radiation transport in these 
simulations. We will combine radiation transport hydro- 
dynamics simulations with non-equilibrium treatment of 
ionization balance in the metal species to study the effect 
of local, non-uniform radiation on WHIM observables. 
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